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Abstract 

In this paper, we present a probabilistic numerical algorithm combining dynamic program- 
ming, Monte Carlo simulations and local basis regressions to solve non-stationary optimal mul- 
tiple switching problems in infinite horizon. We provide the rate of convergence of the method 
in terms of the time step used to discretize the problem, of the size of the local hypercubes 
involved in the regressions, and of the truncating time horizon. To make the method viable for 
problems in high dimension and long time horizon, we extend a memory reduction method to 
the general Euler scheme, so that, when performing the numerical resolution, the storage of the 
Monte Carlo simulation paths is not needed. Then, we apply this algorithm to a model of opti- 
mal investment in power plants. This model takes into account electricity demand, cointegrated 
fuel prices, carbon price and random outages of power plants. It computes the optimal level of 
investment in each generation technology, considered as a whole, w.r.t. the electricity spot price. 
This electricity price is itself built according to a new extended structural model. In particular, 
it is a function of several factors, among which the installed capacities. The evolution of the 
optimal generation mix is illustrated on a realistic numerical problem in dimension eight, i.e. 
with two different technologies and six random factors. 



1 Introduction 

This paper presents a probabilistic numerical method for multiple switching problem with an ap- 
plication to a new stylized long-term investment model for electricity generation. Since electricity 
cannot be stored and building new plants takes several years, investment in new capacities must be 
decided a long time in advance if a country wishes to be able to satisfy its demand 1 . Before the 
worldwide liberalization of the electricity sector, electric utilities were monopolies whose objective 
was to plan the construction of power plants in order to satisfy demand at the minimum cost under a 
given constraint on the loss of load probability or on the level of energy non-served. This investment 
process was called generation expansion planning (GEP). Its output was mostly a given set of power 
plants to build for the next ten or twenty years (see [36] for a comprehensive description of the GEP 
methodology and related difficulties). Despite thirty years of liberalization of the electricity sector, 
of the recognition that GEP methods were inadequate within a market context ([34, 24]) and of an 
important set of alternative methods (see [28] and [22] for recent surveys on generation investment 
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models and softwares), power utilities still heavily rely on GEP methods (see [3]). However, real op- 
tion methods, which should have been the natural alternative valuation method for firms converted 
to a value maximizing objective, did not emerge as the method of choice. Despite the important 
body of literature that followed [45] 's seminal paper, [23] 's monography and implementations for 
electricity generation investment (see for instance [11]), real options still remain a marginal way of 
assessing investment decisions both in the electric sector and in the industry in general (see for in- 
stance, amongst the recurrent surveys on capital budgeting methods, [8]). Nevertheless, as shown in 
[44], firms tend to reproduce with heuristic constraints (such as hurdle rate or profitability index) 
the decision criteria given by real option methodology. 

The main reason for this situation lies in the considerable mathematical difficulties involved in the 
conception of a tractable yet realistic real option model for electricity generation. This difficulty 
reflects in the literature where the main trend consists in designing a small dimensional (1 or 2) 
real option model to assess investment behaviour with respect to some specific variables (see for 
instance [1, 9] for models in dimension 2 analysing the effects of uncertainty and time to build). 
It is still possible to find investment model in dimension 3 based on dynamic programming which 
are numerically tractable (see for instance [46, 11] and Section 5 for comments). But, in higher 
dimension, because of the curse of dimensionality, investment models mainly rely on decision trees to 
represent random factors (see [2] for a recent typical implementation of this approach) . The resulting 
tractability is however obtained at the expense of a crude simplification of the statistical properties 
of the factors. 

Our approach in the present paper takes advantage of the considerable progress made in the last ten 
years by numerical methods for high-dimensional American options valuation problems to propose a 
probabilistic way to look at future electricity generation mixes. For an up-to-date state of the art on 
this subject, the reader is referred to the recent book [15]. 

In this paper, we first adapt the resolution of American option problems by Monte-Carlo methods 
([43, 57]) to the more general class of optimal switching problems. The crucial choice of regression 
basis is done here in the light of the work of [13], so as to obtain a stable algorithm suited to high- 
dimensional problems, aiming at the best possible numerical complexity. The memory complexity, 
often acknowledged as the major drawback of such a Monte Carlo approach (see [16]), is drastically 
slashed by generalizing the memory reduction method from [18, 19, 20] to any stochastic differential 
equation. We provide a rigorous and comprehensive analysis of the rate of convergence of our algo- 
rithm, taking advantage of the works of, most notably, [12], [55] and [29]. Note that such features as 
infinite horizon and non-stationarity are encompassed here. Finally, we build a long-term structural 
model for the spot price of electricity, extending the work of [5] and [4] in several directions (cointe- 
grated fuels and CO2 prices, stochastic availability rate of production capacities, etc.). This model 
is itself incorporated into an optimal control problem corresponding to the search for the optimal 
investments in electricity generation. The resolution of this problem using our algorithm is illustrated 
on a simple numerical example with two different technologies, leading to an eight-dimensional prob- 
lem (demand, CO2 price, and, for each technology, fuel price, random outages and the controlled 
installed capacity) . The time evolution of the distribution of power prices and of the generation mix 
is illustrated on a forty-year time horizon. 

To sum up, the contribution of the paper is threefold. Firstly, it provides, for a suitably chosen regres- 
sion basis, a comprehensive analysis of convergence of a regression-based Monte-Carlo algorithm for 
a class of infinite horizon optimal multiple switching problems, large enough to handle realistic short 
term profit functions and investment cost structures with possible seasonality patterns. Secondly, it 
adapts and generalizes a memory reduction method in order to slash the amount of memory required 
by the Monte Carlo algorithm. Thirdly, a new stylised investment model for electricity generation is 
proposed, taking into account electricity demand, cointegrated fuel prices, carbon price and random 
outages of power plants, used as building blocks of a new structural model for the electricity spot 
price. A numerical resolution of this investment problem with our algorithm is illustrated on a spe- 
cific example, providing, among many other outputs, an electricity spot price dynamics consistent 
with the investment decision process in power generation. 

The outline of the paper is the following. Section 2 describes the class of optimal switching problems 
studied here, including the detailed list of assumptions considered. Section 3 describes the resolution 
algorithm and analyzes its rate of convergence, in terms of the discretization step, of the size of the 
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local hypercubes from the regression basis, and of the truncating time horizon. Section 4 details 
the computational complexity of the algorithm, as well as its memory complexity, along with the 
construction of the memory reduction method. Finally, Section 5 introduces the extended structural 
model of power spot price, the investment problem, as well as an illustrated numerical resolution. 
Section 6 concludes the paper. 



Notation 

Here are some notation that will be used throughout the paper: 

• The notation 1 {.} stands for the indicator function. 

• Throughout the paper, C > denotes a generic constant whose value may differ from line to line, 
but which does not depend on any parameter of our scheme. 

• For any stochastic process X = (X s ) s>0 taking values in a given set X, and any (t, x) e R + x X, 
we denote as X l - X = (Xl' x ) s>t the stochastic process with the same dynamics as X, but starting 
from x at time t: X\ ,x = x. 

• For any (a, b) e R x M, a A b := min (a, b) and aV(>:= max (a, b). 

• Vp > 1, the norms ||.|| p and denote respectively the p— norm and the L p - norm: \fx e E™ 
and any E- valued random variable X such that E [\X\ P ] < oo: 

H P = (Er=iN p )* , \\x\\ Lp =E[\x\r* 

We recall that Vp > 1, Vx e E n , ||a;|| < Hx^ < \\ x \\ 



2 Optimal switching problem 
2.1 Formulation 

Fix a filtered probability space (ft, J 7 , F = (J r t) t >o ; i where F satisfies the usual conditions of right- 
continuity and P-completeness. We consider the following general class of (non-stationary) optimal 
switching problems: 



v (t, x, i) = sup E 



[ X f(s,Xl<*,I?)ds-Y^ fc(r n ,Cn) 

* r n >t 



(2.1) 



where: 



X l - X — {X t g X ) s>t is an Revalued, F-adapted markovian diffusion starting from X t = x e R d , with 
generator C. 

I a = (I™) s>0 is a cad-lag, R d -valued, F-adaptcd picccwise constant process. It is controlled 
by a strategy a, described below. We suppose it can only take values into a fixed finite set 

I q = {ii, «2, . . . , i q }, q € N* with i a = ^£ l^j, which means that equation (2.1) corresponds to 

an optimal switching problem. 

An impulse control strategy a corresponds to a sequence (t„, tn)„ e N of increasing stopping times 
t„ > 0, and T Tn -measurable random variables t n valued in I q . Using this sequence, I a = (/°) g>0 
is defined as follows: 

If = L 1 {0 < S < T } + ^ tnl {Tn<S< T„+l} € Iq 
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Alternatively, a can be described by the sequence (r„, £„)„ gN , where £„ := t„ — i n _i (and Co : = 0). 
Using this alternative sequence, J a can be written as follows: 

J s = k> + E ^ G I, 

• .A is the set of admissible strategies: a strategy a belongs to A if r„ — > +00 a.s. as n — > 00. 

• For any (t,i) G M + x I 9 , the set _A t ,i C A is defined as the subset of admissible strategies a such 
that If = i. 

• / and k are R- valued measurable functions. 
2.2 Assumptions 

We complete the above formulation with the following relevant assumptions. 

Assumption 1. [Diffusion] The Mr -valued uncontrolled process X is a diffusion process, governed 
by the dynamics 

dX s = b(s,X s )ds + a(s,X s )dW s (2.2) 
x Q = x e R d 

where W is a d- dimensional Brownian motion, and b and a are respectively M. d -valued and ~R dxd - 
valued functions. 

Assumption 2. [Lipschitz] The functions b : R + x R d -4 R d and a : R + x R d -4 R dxd are Lipschitz- 
continuous (uniformly in t) with linear growth: 3Cb, C a > s.t. Vi G R+, V (x,x') G : 

\b(t,x)-b(t,x')\ < C b \x-x'\ 

\b(t,x)\ < C b (l+\x\) 

\a(t,x) - a(t,x')\ < C a \x-x'\ 

W(t,x)\ < C a (l + \x\) 

Remark 2.1. Assumption 2 is sufficient to prove the existence and uniqueness of a strong solution to 
the SDE (2.2) (see for instance Theorem 4.5.3 in [37]). 

Remark 2.2. Under Assumption 2, there exist, for every p > 1, positive constants C p and p p such 
that V.s > t > and Vx G R d : 

E[|Xff] <C p (l + \x\ p )e*p( Pp ( S ~t)) (2.3) 

(use Burkholder-Davis-Gundy inequality and Gronwall's Lemma, see for instance [37] Theorem 4.5.4 
for the even power case). 

Assumption 3. [Lipschitz&iDiscount] The functions f and k decrease exponentially in time: 3p > 
s.t. V(t,x,i,j) G R+ x R d x (I,) 2 : 

f(t,x,i) = e-*f(t,x,i) 
k(t,j-i) = e-*k(t,j-i) 

where the functions f and k are Lipschitz continuous with linear growth: 
3C f ,C k > s.t. V{(t,x,i,j),(t,x',i',j')} G [R + xR d x (I,) 2 }': 

\f(t,x,i)-ftf,x',i')\ < C f (\t-l/\ + \x-x'\ + \i-i'\) 
\f(t,x,i)\ < C/(l + |a;|) 
\k^ J - l) ^k (t \j'~ i ')\ < Ck (\t-t'\ + \(j-i)-(j'-i')\) 

Moreover, we assume in the following that p > pi where p\ is defined in equation (2.3). 
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Assumption 4. [Fixed costs] The cost function k : R+ x R d — > M + is smc/i iftai: 
. V< € M +) fc(t,0) = 0. 

. 3 K > o s.t. vt g R+, v(»,j) g (i g ) 2 , ^ i} =► {fc(t, j -*)>»}. 

• (triangular inequality) Vi G M +; V («, j, k) G (1 9 ) 3 tofift i ^ j and j 7^ fc; 

fe(t, k-i) < k(t, j -i) + k(t, k-j) . 

Remark 2.3. The economic interpretations of Assumption 4 are the following: 

1. There is no cost for not switching, but any switch incurs at least a positive fixed cost. 

2. At any given date, it is always cheaper to switch directly from i to k than to switch first from i 
to j and then from j to k. 

Remark 2.4. Under those standard assumptions, the value function v is well-defined and finite. 
Indeed, using equation (2.3), V (to, t, x, i) G R+ x R + x R d x R d ' with t < t and Ma G -4* ,i : 



E 



/OO "I /»OC 

|/( S ,X*«^,/«)|rf S < C> y e-" s (l + E[|A*^|])d S 

< C/ (e- pt + {l + \x\) e- ps e Pl[s - to) ds^j 

< C f (l + \x\)e- pt - plt ° (2.4) 

where p := p — pi > (Assumption 3). In particular, the costs being positive (Assumption 4), and 
recalling (2.1), it holds that: 

v(t,x,i) < C f (l+ \x\)e~ pt (2.5) 
2.3 Outline of the solution 

From a theoretical point of view, the value functions v% :— v (., i G I 9 from equation (2.1) are 
known to satisfy (under suitable conditions on := f(.,.,i) and k, see for instance [53] in 

a much more general setting) the following Hamilton- Jacobi-Bellman Quasi- Variational Inequalities 
(HJBQVI): V (t, x, i) eR+xR d x I 



dvj 

at i&p,jiti 



min <J — (t, as) - £i>i (*,#) - /j (t,a;) , U; (f, a;) - max t (vj (t, x) - k (t,j - i)) \ = Q (2.6) 



together with suitable limit condition. 

Alternatively, the process v(t,X t ,i), t > can be characterized as the solution of a particular 
Reflected Backward Stochastic Differential Equation ([33, 25]). 

Moreover, the value function (2.1) satisfies the well-known dynamic programming principle, i.e., for 
any stopping time r > t: 



v (t, x, i) — sup E 



[ T f(s,Xi'*,I?)ds- ]T fc(r n ,Cn)+«(r,JC*^,Z?) 



t<T,,<T 



(2.7) 



From a practical point of view, apart from a few simple examples in low-dimension, finding directly the 
solution of the HJBQVI (2.6) is usually infeasible, and the numerical PDE tools become cumbersome 
and inefficient in the multi-dimensional setting. Instead, probabilistic methods based on (2.7), in the 
spirit of [16], are usually more practical and versatile. 

Indeed, as the diffusion X is not controlled, this optimal switching problem can be seen as an extended 
American option problem. This suggests that, up to some adjustments, the probabilistic numerical 
tools developed in this context (see [13] for instance) may be adapted to solve (2.1). 
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To be more specific, consider a variant v of (2.1) such that the switching decisions can only take place 
on a finite time grid II = {t = < t x < . . . < t m = T} for a fixed T > 0. Then Mi e l q , Mx e R d , 
and Mtk S II, the dynamic programming principle (2.7) becomes: 



where: 



Ei{T,x) :=E 
Ei(t k ,x) :=E 



max <^ Ei (t k , x) , max {vj (t k ,x) -k(t k ,j — i)} 



fi(s,Xt 



dt- 



k = m — 1, . . . ,0 



(2.8) 

(2.9) 
(2.10) 



and where the notation X t,x := (Xl' x ) s>t refers to the process X conditioned on the initial value 
X t = x. 

If, moreover, the cost function k is such that at most one switch can occur on a given date t k 
(triangular condition), then equation (2.8) can be simplified into: 



Vi (tk,x) = max [Ej (t k ,x) - k(t k ,j - i) Uj^a} 

J Gig 



(2.11) 



which is explicit in the sense that v, (t k , .) directly depends on v, (t k +i, .). 



In practice, apart from the potential approximation of the stochastic process X and of the final values 
(2.9), the difficulty lies in the efficient computation of the conditional expectations (2.10). 

In the American option literature, various approaches have been developed to solve (2.11) efficiently 
Notable examples are the least-squares' approach ([43, 57]), the quantization approach and the Malli- 
avin calculus based formulation (see [13] for a thorough comparison and improvements of these 
techniques). In the spirit of [17], one may also consider non-parametric regression (see [38] and 
[56]) combined with speeding up techniques like Kd-trees ([32, 40]) or the Fast Gauss Transform 
([61, 47, 50, 54, 51]) in the case of kernel regression. 

Here, we intend to solve (2.1) on numerical applications which bears the particularity of handling 
stochastic processes in high dimension (dim (X) = d ^> 3, with however dim (I) = d' ps 3, see Section 
5). For such problems, the most adequate technique so far seems to be the local regression method 
developed in [13]. We are thus going to make use of this specific method to solve (2.11) in practice. 

In the following, we provide a detailed analysis of the above suggested computational method. 



3 Numerical approximation and convergence analysis 

This section is devoted to the precise description of the resolution of (2.1), along the lines of the 
discussions from Subsection 2.3. Moreover, the convergence rate of the proposed algorithm will be 
precisely assessed. 



3.1 Approximations 

Recall equation (2.1) defining the value function v (t,x,i) : 



v (t, x, i) — sup E 



/>oo 

/ f(s,X!?*,I?)d8- J2 k(r n ,Cn) 
Jt _ ^ 



We are going to consider the following sequence of approximations: 
• [Finite time horizon] The time horizon will be truncated to a finite horizon T. 



(3.1) 
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• [Time discretization] The continuous state process X and investment process / will be discretized 
with a time step h. 

• [Space localization] The R d - valued process X will be projected into a bounded domain V e , pa- 
rameterized by e. 

• [Conditional expectation approximation] The conditional expectation involved in the dynamic 
programming equation will be replaced by an empirical least-squares regression, computed on a 
bundle of M Monte Carlo trajectories, on a finite basis of local hypercubes with edges of size 8. 

The rate of convergence of the algorithm will then be provided, as a function of these five numerical 
parameters: T, h, e, M and S. 



3.1.1 Finite time horizon 

The first step is to reduce the set of strategies to a finite horizon: 



VT(t,x 7 i) — sup E 

a€AT, 



9f(T,x,i) :-- 



E 



f f(s,Xl'*,I«)ds- k(T n Xn)+9f(T,X^,I^) 

t<T n <T 

/•OO 

J f{s,Xj'*,i)ds 



(3.2) 
(3.3) 



where < t < T < +oo, and Aj i C A t .i is the subset of strategies without switches strictly after 
time T. Hence the final value gj corresponds to the remaining gain after T. 

Alternatively, one may choose, for convenience, another final value g instead of as long as it is 
Lipschitz-continuous and satisfies a suitable condition (cf. equation (3.21)). The set of such functions 
will be denoted as 6 g/ . The difference between the two value functions is quantified in Proposition 
3.1. 

This freedom on the final values will be used in practice to avoid a computation on an infinite interval 
[T, oo [ as in the definition of gf. 

From now on, we choose and fix one such g £ Q gf . 



3.1.2 Time discretization 

Then, we discrctizc the time segment [0, T]. Introduce a time grid II = {t = < t\ < . . . < tjv = T} 
with constant mesh h. Consider the following approximation: 



vn (i, x, i) = sup E 



f f{s,Xl' x ,If)ds- J2 k(T n ,(n)+9{T,X¥,I%) 



t<r„<T 



(3.4) 



where A] l i C Aj t is the subset of strategies such that switches can only occur at dates r„ £ II fl [t, Tj. 

Now, with a slight abuse of notation, we can safely switch from the notation a = (r n , Cn) n >o to the 
notation a = (T n ,i n ) n>0 (remember Subsection 2.1), replacing the quantity J2t<r n <T ^ ( T n, C«) by 

Et<r„<T fc ( T nJr n -hJ?J or b Y Et< T „<T k ( r « > L n-i , <<„) , where k(t,i,j) = k(t,j - i) . The error 
between w T and wn is quantified in Proposition 3.2. 

Next we also approximate the stochastic process X by its Euler scheme X = (A > s ) 0<s<T , with 
dynamics: 

dX s = fe(7r( S ),A > , (s) )d S + ( r(7r( S ),A > T(s) )dVK s , 0<s<T (3.5) 
X = x £ R d 
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where Vs € [0, T], 7r (s) := max {i £ II; f < s}. The new value function reads: 



vyi (t, x, i) = sup E 



f T f(ir(8),X t a '',If)ds- Yl fc(r„,tn-i,tn)+fl(r s ^*,i?) 



t<T n <T 

The error between i>n and «n is computed in Proposition 3.3. 



(3.6) 



3.1.3 Space localization 

Fix e > 0. Vi € [0, T] , let Z>f be a bounded convex domain of M. d . In particular there exists C (t, e) > 



such that Vr £ Z>£, |a;| < C(t,e). Let 
chosen such that Vs <E [0,i], 

E[|X S -P|(X S )|] <e 



denote the projection on T>\. This domain is 



(3.7) 



Denote this projection as X\ 



XI := VI (X t ) 



In other words, XI is equal to X t most of the time (i.e. when X t € X>f), except when X t is outside 
2?( , in which case Xf corresponds to the projection of X t onto T>f. 

Define vfi as ^he value function vu from equation (3.6) with X replaced by X e . The error between 
those two value functions is computed in Proposition 3.4. 

Example 3.1. To clarify this construction of space localization, we explicit it on the very simple 
example of a d-dimensional standard brownian motion (Wt) t£ r 0T i. In this case, X t = X t — Wt- 

Choose T>\ to be a centered, symmetric hypercube: T>\ — [— C (t, e) , C (t, e)] d for some constant 
C (t,e). Hence, Va; £ R d , VI (x) := —C(t,e) AiVC(t,£) component-wise. With this expressions, 
one can find a C (t,e) such that (3.7) holds. Indeed, Vs € [0,7*]: 

E [\W S - VI {W s )\] < E [\W t - VI {W t )\] = E [\W t -C(t,e)\ l {lWtl>C (t,e)}} = 2 d E [(W, 1 - C(t,e)) 

(3.8) 

where W 1 is a one-dimensional Brownian motion. Hence, finding a value for C (t,e) such that (3.7) 
holds boils down to inverting Bachclicr's option pricing formula in order to get the strike as a function 
of the price of the call option. This is done in [6], see [52], but under the form of a series expansion 
for small moneyness, which is unsuitable for our purpose (because C (i,e) — > oo when e — > 0). Thus, 
we are here only going to look for a simply invertible upper bound for (3.8). Denoting as Af the 
cumulative distribution function of a standard Gaussian random variable, and using the standard 
inequality 1 — Af (x) > 



E 



(W? - K) 



(xy/t - K^j Af' (x) dx 



Vi 

y/2tt 



K l-Af 



K 
Vt 



< 



K 2 



K 2 + t 



2t < . e 2t 



2tt 



Inverting this last upper bound, the inequality (3.7) is satisfied with C (t,e) 



3.1.4 Conditional expectation approximation 

From now on, in order to prevent the notation from becoming too cumbersome and clumsy, we are 
going to drop the e index in the following, i.e. X t will stand for XI, and vn for vfi- 

For the fully discretized problem (3.6), the dynamic programming principle (2.11) becomes: 



Vn(T,x,i) =g(T,x,i) 

v u (*n, x,i) = max {hf(t n ,x,j) - k (t n ,i,j) +E v u (t n+ i, X 1 ^, j^j j ,n=N-l,. 



,0 (3.9) 
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The last step is to approximate the conditional expectation appearing in equation (3.9). As discussed 
in Subsection 2.3, we choose to approximate it using the following regression procedure. Consider 
basis functions (eu (x)) 1<k<K , K G N U {+00}, x S R d . For suitable functions <p : II x R d x I q — >• E, 
define: 



A = A," (ip) := are min E 
1 \eR K 



-Z) A * efc (^*») 



(3.10) 



As truncating the approximated conditional expectations is a necessity in theory as well as in practice 



(see [12, 30, 55]), suppose that there exist known bounds T tn ' x (y>) and T "' (ip) on E 



r^^)<E^ (t n+ll I t ^,i) 



Then, Vi G I q the quantity E 



is approximated by: 



95 [t n+1 , Xfc*,i 
(3.11) 



E 



AT 



fc=i 



which is used to define the next approximation vn of the value function: 
v n (T,x,i) = g(T,x,i) 



(3.12) 



v n (t n ,x,i) = max j/i/ (f n , x,j) - k (t n ,i,j) + E w n (tn+i, x t2+i>j) } > = N - 1, . . .QD13) 



Interesting discussions on the choice of function basis can be found in [13]. In particular they advocate 
bases of local polynomials, which is numerically efficient and well-suited to tackle large-dimensional 
problems (see Subsection 4.1). However, for the sake of simplicity, we will restrict our study in this 
section to a basis of indicator functions on local hypercubes (as in [55] and [30]) (which is the simplest 
example of local polynomials) as defined below: 

For every t n G II, consider a partition of the domain X>f into hypercubes \B^ J 



k=l, 



and Bl 



n B l = 



Vz ^ j. It may be deterministic, or J^-measurable. We 



only assume that there exists (6, S) € with 5 < 5 such that the lengths of the edges of the 
hypercubes, in each dimension, belong to [5, S] (in particular, the volume of each hypercube B% 



belongs to 



5 d ,5 d 



) . This liberty over the definition of the partition enables to encompass the kind 

of adaptative partition described in [13]. Then, the basis functions considered here are defined by 
e k tn (x) :=l{i£ B* n }, x G R d , 1 < k < K e . 

With this choice of function basis, the error between vjj and vu is computed in Proposition 3.5. 

Finally, let (^t™))^™^"^ be a finite sample of size M of paths of the process X. The final step is to 
replace the regression (3.10) by a regression on this sample: 



M 



Then Vi G I„ the quantity E 



(p (t n+ i,X*"'*,i 



<p(t n+ i,X™ +1 ,i 




K 

^ Afcefe (XI 

k=l 



is approximated by: 



(3.14) 



E 



K 



E^^V^W (x) AT tn ' x (ip) 



fc=l 



(3.15) 



leading to the final, computable approximation Dn of the value function: 
vu(T,x,i) = g(T,x,i) 

v n (t n ,x,i) = maxlhf(t n7 xJ)-k(t n ,iJ) + Ev n {tn+i,X t t ^ i ,j) \ , n = N - 1, . . .Q216) 



The error between vjj and vjj with the same choice of function basis is given in Proposition 3.6. This 
proposition will make use of the following quantity: 

p(t n ,S,e) := min min P (X t £ B?) (3.17) 

tenn[o,t„] B*cz>f 

which is strictly positive, as the domains T)\, t € [0, T] are bounded. More precisely, only lower 
bounds of these quantities will be required. 

Example 3.2. Carrying on with Example 3.1 of a d-dimensional Brownian motion, we explicit a 
lower bound for p(t n ,S,s) in this simple case. First, P (W t 6 B* ) = f Bk fw t (x) dx where fw t is 

the density of W t . As Vfc, Bj? n C V% n , where in this example V\ = [-C (t, e) , C (t, e)] d with C (t, e) = 
Jtm(-%), it holds that V.t e 2?f, /w t (*) > (/W? = pfp Hence P (W t € B t fc ) > 

Tj^Vol(-B^) > ^|ya-5 d . As a conclusion, p(t n ,8, e) > . d S . Remark however that this lower 
bound is very crude, and that it can be very far below p (t n , S, e) for large S. 



Combining all these results, we obtain a rate of convergence of z)n towards v: 



Theorem 3.1. \/p > 1 , 3C P > such that: 



max|u (t ,x ,i) - v u (t ,x ,i)\ 



< C p { (1 + |x |) e-P T + (1 + |z |) s Vft + e + ^ + 



l + C(T,e) 



, i + g(r, £ ) 

/i>/Mp (T, 5, e) 1_ ^ ( T > ^ £ ) 



In particular, vu (0,Xo,i) — >l p v (0,Xo,i) uniformly in i e I, when T — > oo, h — > 0, e — > 0, S — > 
and Af->oo with f -> 0, 1+C(T ' £ 1 ) _^ L -4 and 

h%/Mp(T,<5,e) pv2 



l+C(T,e) 
hMp(T,8,e) 



0. 



Remark 3.1. If the cost function fc (recall Assumption 3) were to depend on x, then, under a usual 
Lipschitz condition on k (similar to that of /), Theorem 3.1 would still hold, replacing only the term 

(1 + |sco|)^ V^by (l + \xq\^\ y h log (^) (recalling Remark 3.4). 

Remark 3.2. The adaptative local basis can be such that each hypercube contains approximately 
the same number of Monte Carlo trajectories (see [13]). This means that p ^ s ^ ~ b where b is the 
number of functions in the regression basis. With this remark in mind, the leading error term in 
Theorem 3.1 behaves like for p = 2. This is close to the corresponding statistical error term in 

[42] 61 ° g ff f) ) in the context of BSDEs. The advantage of their approach is that they can handle 
any (orthonormal) regression basis, while our approach (in the context of optimal switching) provides 
a bound on the L p error for every p > 1. 

Example 3.3. In the case of a d-dimensional Brownian motion, the rate of convergence of Theorem 
3.1 can be explicited further, using the upper bound on C(T,e) from Example 3.1 and the lower 
bound on p (T, 5, e) from Example 3.2. Moreover, one can express the rate of convergence as a function 
of only one parameter, choosing the five numerical parameters T, h, e, S and M accordingly. For 
instance, assuming 5 = 6, and minimizing over 6, h, e and T, one can get a convergence rate upper 
bounded by C p (1 + |a:|) 5 z by choosing M ~ z _ 5[ 6 ( d + 1 )l . This is admittedly highly demanding in 
terms of sample size M, but remember that this expression suffers from the crude lower bound on 
p (T, S, s) we chose previously. 



3.2 Convergence analysis 



From now on, we suppose that all the assumptions from Subsection 2.2 are in force. 
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3.2.1 Finite time horizon 

Lemma 3.1. There exists C > such thatV(t,x,i) E R+ x R d x R d ' : 

< v (t, x, i) - v T (t, x,i)<C(l + \x\) e-^vr-^t 

Proof. First, we introduce the following notations: 



(3.18) 
(3.19) 

(3.20) 



H(t,T,x,a) := f f {s,X\<* ds - £ fc(r„,C„) 

* t<r„<T 

J(t,T,x,a) := E [if (i, T, a, a)] 
for any admissible strategy a e ^4 ti j. In particular: 

v(t,x,i)~ sup J (t, +oo, x, a) , «t (i, x,i) = sup J(i, +oo,x,a) 

Fix (t,x,i) e R+ x R d x K d '. Using equation (3.20): 

VT(t 1 x,i)= sup J (t, co, a;, a) < sup J (t,oo,x,a) = v (t,x,i) 

aeAf i ot£A t ,i 

which provides the first inequality. Consider now the second inequality. Choose e > 0. From the 
definition of v (equation (3.1)) there exists a strategy a £ <G At.% such that: 

v (t, x, i) — e < J (t, oo, x, a e ) < v (t, x, i) 

Define the truncated strategy a e T £ Af A such that Vs e [t,T], lf T = if and Vs > T, Jf = if. In 
order not to mix up the variables r„ and ( n from different strategies, we add the name of the strategy 
in index when needed. Then: 

H (t, oo, x, a £ ) — H (t, oo, x, a?) 

= f( s ,xr,if)ds- k(r:\c)^ 

rOO 

/ f( s ,xr,if)d S - k(r:Kcn T ) 

J™ f(s,X^,lf)ds- g fc(rf,Cf)| 

- J r vT / ( s , if) d S + r f (a, xr,if T ) ds- yi k {< > c e ) I 

= r f(s,x^,if)ds- r f(s,x^,if T )ds~ k (<^n) 

< f° / (s, JtfMf ) ds- r f (s, X*;*, lf T ) ds 

as k(s, 0) = and fc > (Assumption 4). Hence, using Jensen's inequality and equation (2.4), 
3C > such that 



< 



| J (i, oo, x, a e ) — J (t, oo, x, a^)| < E (i, oo, x, a e ) — H (i, oo, x, af*) 



< E 



/"°° / (s, X^, if) da] +E [ r f (s, X^,lf T ) 



< C(l + \x\)e 



-ptVT-pxt 
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Finally, given that v(t, x,i) < e + J (t, oo, x, a £ ) and (t, x, i) > J (t, oo, x, a T ) , the following 
holds: 

v (t, x, i) — vt (t, x, i) < e + J (t, oo, x, a 6 ) — J (t, oo, x, a T ) 
< e + C(l + |a;|)e-^ vT -' ,lt . 

Since this is true for any e > 0, and that C, p and p\ do not depend on e, the proposition is 
proved. □ 

Now, we focus on the final boundary gf. For the time being, denote the value function (3.2) as v T f 
to emphasize the dependence of v on the terminal condition. As a consequence of equation (2.4), 
V(x,i) e R d x I,: 

| fl/ (T,a:,i)|<C(l + |x|) e -^ (3.21) 



Hence, define the class Q gj of Lipschitz functions from 

V (T, x, x', i)eI+xR d xR' l x I q : 



x I„ into K such that \/g € 6 



9f- 



\g(T,x,i)-g(T,x',i)\ < C g e~r T \x-x'\ 
\g(T,x,i)\ < C g e-" T (l + \x\) 



(3.22) 
(3.23) 



for some C g > 0. Obviously <?/ G g/ . Now, for any g S 6 g/ , denote as the value function defined 
as in equation (3.2) with g instead of gf. We are going to show that the precise approximation error 
due to the choice of final value g does not matter much as long as g is chosen in this class ff/ . 



Lemma 3.2. There exists C > such that V (t, x, i) £ R + x R d x l q : 



v T f (i, x, i) - v 9 T (i, x,i)\ < C (1 + \x\) e 



-ptvT-pit 



Proof. Fix (t,x,i) 6 M + x R d x I,. To shorten the proof, we assume that (resp. v T ) admits an 
optimal strategy a*j € Aj t (resp. a* £ Aj ^ (this assumption can then be relaxed using e-optimal 
strategies as in the proof of Proposition 3.1) 1 . Therefore, recalling the notations H (equation (3.18)) 
and J (equation (3.19)) introduced in the proof of Lemma 3.1: 



xSp (t, x, i) — v T (t, x, i) = J (t,T, x, off) + E gj (t, 



— J(t,T,x, a*) — E <?(T,X^V£ 



j(t,T,x,a* f ) +E g(T,X^ x ,I^ — J (t, T, x, a*) — E g{r,X^ x J^ 



g f (T,X^ X ,I^~ 9 (T,X T ^,I^ 

g f (T,x¥,4 } )-g(T,x¥,I%' 
< C(l + E[\X^ x \])e- pT < C(l + e" 



< E 



ptvT-ptt 



Symmetrically, the same inequality holds for v T (t, x, i) — v T (t, x, i), ending the proof. 



□ 



Proposition 3.1. There exists C > independent of T such that V(t,x,i) £ R+ x R d x I q and 
Vg £ e 



9/ • 



\v (t, x, i) - 4 (i, x, i)\ < C (1 + \x\) e -^ VT -^ 4 
Proof. Combine Lemmas 3.1 and 3.2. 



□ 



From now on, we choose and keep one final value function g £ Q g , , and remove the index g from the 
notation of v and its subsequent approximations. 

^^Note that under the assumptions from Subsection 2.2, one may use Theorem 3.1 from [35] to get the existence of 

2" 

a unique optimal strategy a* for the value function (3.2), satisfying E X/0<r a * <T ^ ( T ™ ' ^" ' x 
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3.2.2 Time Discretization 

Proposition 3.2. There exists a positive constant C such that for any (t,x,i) € II X M. d x I q : 

\v T (t,x,i) -v n (t,x,i)\ < Ce- pt (l+ \x\^ hh (3.24) 

Proof. Under the assumptions from Subsection 2.2, one can apply Theorem 3.1 in [29] to prove (3.24), 
noticing that the cost function k does not depend on the state variable x. 

Use the discounting factor in the definition of / to factor the e~ pt term and to get that C does not 
depend on T. □ 

Remark 3.3. Another alternative to get this rate of h^ is to work with the reflected BSDE represen- 
tation of Vji, as in [16] (adapting [12]) or [21]. 

Remark 3.4. Were the cost function k to depend on the state variable, the upper bound in Proposition 
3.2 would only be Ce~ pt (l + \x\^ (ft, log (^)) 5 , as stated in [29] (making use of results from [27]). 

Proposition 3.3. There exists C > such that for any (t, x, i) S II x W 1 x I q : 

\vn (t,x,i) - v u {t,x,i)\ < Ce~ pt h^ 

Proof. T and g being fixed, we can define, in the spirit of equations (3.18) and (3.19), the following 
quantities: 



H (i, x, a) 

J (t, x, a) 
H (i, x, a) 

J (t, x, a) 



f f(s,Xl' x ,I?)ds- MWn-i,^) + <7(*VT,X t ^ T ,J t V) (3.25) 

t<T„<T 

E [H (t, x, a)] (3.26) 
I f(n(s),Xi< x ,I?)ds- ^(T„, l „-i,^)+ 5 (tVT,X t t ; :c T ,/« T )(3.27) 

* t<r„<T 

E[H(t,x,a)] (3.28) 



for any admissible strategy a € A^i- For these discretized problems, the existence of optimal controls 
a* and a* is granted. Hence: 



i>n (t,x,i) - v n (t,x,i) 



= J (t,x,a*) — J (t,x,a*) 

— J (t, x, a*) — J (t, x, a*) + { J (t, x, a*) — J (t, x, a*)} 
< J (t, x, a*) — J (t, x, a*) 



e- ps E 



ds 



f[s,Xl\I« j-f{n(s),X^,C 
(T,X^ x ,I^-g(T,X^,I^ 

r T 

< C f J e~ ps E [\X l s ' x - X\' x \] ds + C g e- pT E [\X^ X - X^ x \] 



-E 



< Ce~ pt E 



sup \Xt' x - 




t<s<T 





< Ce~ pt hh 



using the strong convergence speed of the Euler scheme on [i,T]. Symmetrically, the same inequality 
holds for vn (t, x, i) — vjj (t, x, i), ending the proof. □ 



3.2.3 Space localization 

Recall from Subsection 3.1.3 the definition of the bounded domain If, t £ [0,T]. 
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Proposition 3.4. Ve > 0, there exists C > such that for any (x,i) £ R d x I, 

\vu {0,x,i) - vfj (0,x,z)| < Ce 



Proof. Recall the definitions of H(t,x,a) (equation (3.27)) and J(t,x,a) (equation (3.28)), and 
define the quantities H e (t, x, a) and J e (t, x, a) like H (t, x, a) and J (t, x, a) but with X replaced by 
X e . Then, for every (t, x, i) G II x R d x l q and a G Af s : 

J(t,x,a) = J s (t,x,a)+ J T E[f(ir(s),Xi<*,I?)-f(TT( S ),XI^,I?)]d S 
+ E[g(T,X¥,I%)-g(T,X s /< x ,I%)] 

But: 

£ E [/ (tt (s) , Xt>*,I?) - f (tt (s) , ds + E [g (T, - g (T, J?)] 

< C/ e-""E [|X*' X - rfs + C g e-" T E - X £ /' x \] 

It follows that: 

\v u (t,x,i) -v £ n (t,x,i)\ < C f J e- ps E [\X^ X -X e /' x \] ds + C g e- pT E [\X^ X - X^' x \] 

In particular, at t = 0, using equation (3.7), 3C > such that: 

\vn {0,x,i) - «n (0,x,i)| < Ce 



□ 



3.2.4 Conditional expectation approximation 

From now on the domains T>f,t£ [0, T] are fixed once and for all, and, with a slight abuse of notation, 
we will drop e from the subsequent notations. 

We start with preliminary remarks. Recalling Subsection 3.1.4, with this choice of basis, A*™ (ip) 
(equation (3.10)) and A*™ (ip) (equation (3.14)) become: 

~ t , E\(p(t n+1 ,x t ^,i)i\x t e-BMl r . - .. _ 

\'" if) = L { + ' V + ; D K — = E [<p (t n+ i,X tn+1 , i) | X tn G B* ] , 1 < fc < K e 



M 



ez=i ^ [tn + i,xr n+1 ,i) i {x? n G a* } 



Extending these equations, define 

V (w : = p(x t efi, (x)) [<p{t n +i,x tn+1 ,i)\x tn e B tn (x)J (3.29) 



for every (t n ,x, i) G IIx2? n x I 9 , where Vx G , i?t n (x) is the unique hypercube in the partition 
which contains x at time t n , M x n := {m G [1,M] , X™ G i? t „ (x)} and Mf n := #M x n . 

Finally, recalling the approximated conditional expectations (3.12) and (3.15), 
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define for any (t n ,x,j) £ HxT> n x I q and any measurable function p : II x R d x I q — >• R, the following 
quantities: 



$*"' x (^) := E 
:= E 



93 (t„ + i,X^, j 



(t n +i, Xt^'+i 



(3.31) 
(3.32) 
(3.33) 



t n ,X 



where (recalling equation 3.11) T tnm ' x (tp) and r* n '"" (ip) arc lower and upper bounds on 3>*. n,a; (<p): 



Remark 3.5. These definitions are useful to express the dynamic programming equations (3.9), (3.13) 
and (3.16). Indeed, these equations become: 



v n (T,x, 
v n {t n ,x, 

v n {T,x, 
v n (t n ,x, 

v n {T,x, 
v n {t n ,x,i 



= 9 ( T , x, i) 

= m&x{hf{t n ,xJ)-k(t n JJ) + $ t j n ' x (v n )}, n = N -1,...,0 

= 9{T,x,i) 

= max{hf{t n ,x,j)-k{t n ,i,j) + 3>f' x {v n )} , n = N - 1, . . . ,0 

= g{T,x,i) 

= max\hf(t n ,x,j) -k(t n ,i,j) +& n < x (v n ) \ , n = N - 1, . . . ,0 
jei q y. J 



Remark 3.6. For ip — u n , we can easily explicit bounding functions T tn ' x (v n ) and r*"' x (v n ) of 
fy 1 ™'* (vn)- Indeed, using the growth conditions on / and g, the nonnegativity of k and the definition 
of C(T,s) (see Paragraph 3.1.3), there exists C > such that V(t n ,x,j) <G HxV u x I g : 

\v u (t n ,x,j)\ < Ce-"*»(l + C(T,e)) 



|** B,x («n)| < r*»(«n) ~Ce-^(l + C(T,e)) 



(3.34) 
(3.35) 



Moreover, the same is true for p = vn- there exists C > such that V (t n ,x,j) € HxT> n x I q : 

\v u (t n ,x,j)\ < Ce-^(l + C(T,e)) (3.36) 
\¥f' x (vn)\ < r'"(i)n) :=Ce-^(l + C(T,e)) (3.37) 

Finally, we impose the same bound for the definition of vn, i.e. T tn (tin) := r* n (vn)- 

Now we can start the assessment of the regression error. 

Lemma 3.3. Consider a measurable function p : II x R d x I q — > K. Suppose that, for a fixed 
t n+ i e II, it is Lipschitz with constant C n +i, uniformly in j: V (x\,x 2 , j) € K d x R d x I q 



\<p(t n+1 ,x 1 ,j) - <p(t n+1 ,x 2 ,j)\ < C n+ i \xi - x 2 \ 



(3.38) 



TTien &j" x {f) is Lipschitz with constant C n+ \ (1 + Lh), uniformly in j, where L := C& + -f > 0. 
Proof. Choose (f n , j, Xl , x 2 ) e II x I, x M d x K d . Then: 



*5" X1 (v)-$5" ,xa (p) 



3 



E 



< 
< 



p (t n+1 , x\ n n ^ ,i)-ip (t n +i , x\ n n ^ , 

</> {tn+1 , X^ 1 ,j) -<P (t n +l , X^ 2 , j) 
(t n+1 , X^ 1 , j) - (t n+1 , , j) 
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Now, using equations (3.38) and (3.5), and G denoting a d-dimensional standard Gaussian random 
variable, we have 



E 



(y> (t n+ i,x t 



t n ,X\ 

+1 



,j) - P (t n +l,X t 



+ 1 



( Yt n ,X! v t n ,x 2 \' 

[ A t n+1 - A t n+1 ) 

(xi - x 2 + h(b(t n ,x 1 ) - b(t n ,x 2 )) + Vh(<r (t n ,xi) - cr(t n ,x 2 )) G^j 
{(^l - x 2 + h(b(t n ,xi) -b(t n ,x 2 ))) 2 + hE ((a(t n ,xi) -<j(t n ,x 2 )) G) 2 J j 



< Cl +1 {x x - x 2 ) 1 {1 + (2C b + Cl)h + C 2 b h 2 } 



< C 2 n+1 (x 1 -x 2 ) 2 (C b + 



c b + 



+ h\ . 



Thus: 



x, {(p) ^ X2 {v) \< Cn+1 li + l Cb +^\h)\x 1 - x 2 



\t n ,x 2 



Cl 



□ 



Lemma 3.4. Consider again a function ip : II x R d x I„ -)• R such that (3.38) holds for a given 
t n+1 e n. Then, V(x,j) e V £ tn x I,: 



In particular: 



|A*».« ( V ) - (^)l < C n+ iJ(l + Lh) . 



!$*».* ( V ) _ (^)l < C n+1 <5(l + L/i) 



(3.39) 



Proof. Recalling the definitions of B tn (x), of \™' x (tp) (equation (3.29)) and of $f' x (<p) (equation 
(3.31)), simply remark that: 



_ min <J>*"' £ (ip) < $* n,x ((p) < _ max (<p) 



iefl t „(i) 



min (p) < X" ,x {<p) < max $ (<p) . 



x£B tn (a 



ieB t „(i) 



Now, using Lemma 3.3 



A—M-^-^M < _ max &f' x {ip) - min <f>f' x (ip) 

< C n+1 (l + Lh) 

(a 

< C n+1 {l + Lh)5 



< C n +i (1 + L/i) max |a;i — x 2 \ 

(xi,x 2 )e-B t „(a;) 2 



Lemma 3.5. V(i n , zi,x 2 ,i) S n x x I g : 

\vn (t n ,x\,i) - v n {t n ,x 2 ,i)\ < C n \x\ - x 2 \ 

where: 



Cn 
C n 



- ptN C„ 



hCfe~ ptn + C„+i (1 + Lh) , n = N - 1, . . . , 
In particular, 3C > such that Vn = 0, 1, . . . , iV: 

C n < Ce-^e^-^ 



□ 



(3.40) 



(3.41) 



(3.42) 
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Proof. Recall Remark 3.5. We prove the lemma by induction. First, remark that, using hypothesis 
(3.22), it holds for n — N. Now, suppose that it holds for some (n + 1) € [1, ... ,N]. Then, using 
Lemma 3.3: 

vn (t n ,xx,i) 

= max {hf(t n ,xi,j) - k(t n ,i,j) + ^"'^(un)} 

= max {hf(t n ,x 2 ,j) - k(t n ,i,j) + $f' X2 (v n ) + h(f(t n ,x 1 ,j)-f(t n ,x 2 ,j))+ ' Xl («n) ~ ^f' X2 (vu))} 

3 

< max {hf(t n ,x 2 ,j) - k(t n ,i,j) + $f' X2 (v n ) + he~ pt -C ] \x x - x 2 \ + C n+1 (1 + Lh) \ Xl - x 2 \) 

3 

= v n (t„,x 2 ,i)+ (he- ptn C f + C n+ x(l + Lh)) \x x - x 2 \ 

Symmetrically, the same inequality holds for vn(t n ,x 2 ,i) — Vn(t n ,xi,i), yielding equations (3.40) 
and (3.41). Finally, use the discrete version of Gronwall's inequality to obtain equation (3.42) □ 

Proposition 3.5. 3C > s.t. V(t,x,i) enxM^xI, : 

\v n (t, x, i) -v n (t,x,i)\< C-e- pt . 

Proof. For each t n € II, we look for an upper bound E n , independent of x and i, of the quantity 
|% (t n ,x,i) - v u (t n ,x,i}\. First: 

|«n (T,x,i) - v u (T,x,i)\ = \g(T,x,i) - g(T,x,i)\ = 

Hence En = 0. Fix now n £ [0, N — 1]. Using Remark 3.5: 

vu(t ni x,i) = max{hf (t n ,xj) - k(t n ,i,j) +¥ j n ' x (v u )} 



max {hf (t n ,x,j) - k(t n ,i,j) +& n ' x (v n ) 



( % ) - (y n ) 

+^r x (vu) - ^r x (vn)} 

Using Lemmas 3.4 and 3.5, $j n ' x (vn) — (vu) < C n +i6 (1 + Lh) where C n +i is the Lipschitz 
constant of vn at time t n +\ (see Lemma 3.5). Moreover, 

$f' x («n) - («n) < E[%(t n+1 ,4 w) j) -jjn(t n+1 ,I in+1) i)|X t „ eBji)] 
Hence: 

wn (t n ,x,i) < vu (t n ,x,i) + C„+i<5 (1 + i/i) + S n+i 
Symmetrically, the same inequality holds for vu (T, x, i) — vu (t n , x, i), leading to: 

|t>n (t n ,x,i) - vu (t n ,x,i)\ < E n 

Where: E N = 

E n = C n+ i^ (1 + Lh) + -E„+i . 

Consequently, using equation (3.42): 

N x 



S„ = 5(l + L/i) £ C fc <C 



fe=n+l 

where C > does not depend on t n nor T. □ 

The following lemma measures the regression error. It is an extension of Lemma 3.8 in [55] (itself 
inspired by Theorem 5.1 in [12]). 
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Lemma 3.6. Consider a measurable function ip : H x R d x I q — > R. For any p > 1, i/iere ezisis 
C p > sucft thaty (t n ,l,j) e n x [1,M] x I 5 : 



(3.43) 



^ P (X tn e B tn (Xl)) 1 -^ M p e B ^ 

where (p tn € M+ is swc/i £/iat |<p (i n +i, -Xt n+1 , j) | < as. ■ 

Proof. Define the following centered random variables: 

-i 1 M 

(?) : = m £ ^(Wa t : +1 ,j)l{^ G S tn (^J} -E^(t„ +1 ,X- +i ,j)l{X™ e B tB (^J} 

m— 1 

M 



Then: 



< 



a 2r*- O) 



(1) 



<-i+2r*- (^)i 



and: 



£ *n.^„ (l) 



„ 1 



(1) 



1 



'(x tB eB tB (^J) 2 



*7' tn M i 



(1) 



£*n.^n (1) 



< 



'(4^ "2 



< < 



to) 



M 



A 3r*" (?) + 



-tn,Xi 



(1) 



(1) 



M 



< 



< 



P(x t „efl tn (X|J) 

Now, for any p > 1: 



t„.,x 



C' tn to) 



A5r*" (?) + e*""^» (i) r*» (?)ll 



(1) 



< 



'{X tn GB tn (X l t J) "2 



^« (?)-^"'*'" to) 



< 



2 3p- 



-t n ,x' 



er Xtn to) 



(1) 



A5r'"(?)J +{|e*"^(l)|r t "(^)} P |x 

' (1) 



11 I £ 



> 



<(X tn eB tn (XlJ) 2 
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and: 
E 

< 



o3p-2 

' E 



. (X tn e B tn (xlJY 
+ 2 2 p- 1 (r*- (^)) p i 



a 5r*- (ip) 



+ (r*-(^)) p E e 4 "^(l) 



< 



8 p 



(x tn eB tn (xl n )Y 



E 



(1) 



> 



>(x tn eB tn (xl)) p ' 



2p 



a 5r*» (cp) 



+ {r'» m} p e 



using Markov's inequality. Now, the following lemma will provide upper bounds for E 



(3.44) 



(1) 



and E 



Lemma 3.7. For eijery p > 1, there exists C p > smc/i i/iat /or any i.i.d. sample X\, . . . ,Xm of 



\Xi 



ipV2 



< oo, i/ie following holds: 



C p 



, , ■ 1 1 -kpV2 



-valued random variables such that E [Xi] = and E 

1 M 
M E X " 

m—l 

Proof. Using Marcinkiewicz-Zygmund's inequality, there exists C p > such that: 

' M 
E \ X m 



(3.45) 







M 


P- 




E 




E Xm 




< C P E 






m—l 







\m—l 



Multiplying both sides by tt^: 



E 



1 M 
M E X « 



m—l 



Ml 



^ M N 

if E i x ™i 2 



M 



m—l 



(3.46) 



If p > 2, then | > 1 and, using Jensen's inequality: 



Taking expectations on both sides: 



E 



1 M \ 2 M £ M 



M 



m—l 



<E[|X!n 



(3.47) 



Now, if p < 2, then | < 1 and, using Jensen's inequality: 



E 



M N 
1 x-^ , „ ,2 



M 



E^ 



m—l 



< E 



M \ 
1 , — ,2 



M 



m—l 



E 



l*i| 



(3.48) 



Then combine inequalities (3.46), (3.47) and (3.48) and take the power ^ to obtain inequality (3.45). 

□ 
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Now, suppose that 3ip tn £ R + s.t. \ip (t n+ i,X tn+1 ,j)\ < ip tn a.s. . Then, using Lemma 3.7, 3C P > 
such that: 



E 
E 



(1) 



Ml 



1 {Xt n e B tn (X\ n )} - P (X tn £ B tn (X l tn )) 



l \\lfV2 



(*>) 



' c. 



p ' M p 'Ms 



-E [^(t n+1 ,X tn+1 ,j) i{x tn eB tn (x l tn )}] f 2 ' 
where, for the second inequality, the term m — / in the sum was treated separately. Then: 

e [\<p(t n+1 ,x tn+1 ,j) i{x tn e B tn (x\ n )} -E[<p (t n+1 ,x tn+1 ,j)i{x tn e B tn (x l tn )}}f 2 

< (a^E [(^) pv2 i{x tn e B tn (xi)} + E[(^) pv2 i{x tn e B tn (xl n )}]])^ 2 
<V> (^) p P(x tn eB tn (X\ n ))^ 

In a similar manner: 



(3.49) 



(3.50) 



E 



\i{x tn £ B tB (^J} -F(x tn e B tn (xl))f 2 pV2 <»(i t „e5 t „(i t ! J) 



W \ \ P V2 



(3.51) 



(3.52) 



Finally, the combination of inequalities (3.44), (3.49), (3.50), (3.51) and (3.52) proves equation (3.43). 

□ 

We now apply Lemma 3.6 to vu in the following Corollary: 

Corollary 3.1. For every p > 1, there exists C p > s.t. V(t n ,l,j) £ U x [1,M] x I 9 : 



l>*" ,Xt " («n) -l»*"' Xt " (un) 



_ ptn 1 + C(T, £ ) 



1 



Mp(T,5,6)' 



1 

Mp(T, 5,s)' 



Proof. First, recall from equation (3.36) and (3.37) that there exists C > such that for every 

(t n ,j) enxi,: 

Tf(vn) = Ce-*~(l + C(T,e)) 
\vn{t n+ i,X tn+1 ,j)\ < Ce-^(1 + C(T,6)) 

Hence one can apply Lemma 3.6 to vu with these upper bounds. The final step is to recall that the 
minimum probability p(T,5,e) defined in equation (3.17) is a lower bound on P [X tn £ B tn (X l t )) 
for any (t n ,l) £ II x [1,M\. "□ 

Using this result, we can now assess the error between vjj and vu- 
Proposition 3.6. Vp > 1, 3C P > s.t. V Z) £ Id x [1,M] : 

1 + C(T,6) 



sup \v n (t, X\ n ,i)-vn (*, X\ n , i) 

iet n 



< C p e~ pt " 



1 + 



— f 

where I*™ is the set of T tn -measurable random variables taking values in I q . 

Proof. For each t n £ II, we look for an upper bound E n , independent of /, such that: 



sup \v n {t, X\ n ,i) -v n (t, X\ n , i) 



< E„ 
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First: 



sup\v n (T,X^,i)-v„ (T, X l T , \ 



sup\g(T,X l T ,i)-g(T,X l T ,i) 



= 



Hence En — 0. Fix now n £ [0, N — 1]. Recall the dynamic programming equations from Remark 
3.5, and, for every (i,l) £ I*™ x [1,M], introduce j* (resp. j*) the argmax for vn (resp. vn) at point 
X\ n , i.e.: 

v n (t n ,X l tn ,i) = hf(t ni XlJ*)-k(t n ,iJ*)+$f: xL (v n ) 
v n {t n ,X l tn ,i) = hf (t n ,X\J*) - k (t n ,i,j*) + &£ Xln (fin) 



Now: 



v n (t n ,X\ n ,i) = hf (t n ,X l tn ,j*^ - k (t n ,ij*) +^"' Xt " («n) 

/»/ (t n ,X l tn ,j*) -k(t n ,i,j*) +&r XL (vn)} 



^r 1 '" («n) - < f' ; r Xt '* (%) i> + ^ ^r 1 '" m - * 



< 



(%)} 



* (vn) -** B, i " (v n ) 



sup 



$ • n '* t " (%) - *'"' Xt " (Vn) 



Symmetrically: 



«n(t n ,^ B ,j) < v n (t n ,X l tn ,i) + 

jeh 



^"' Xtn (v n ) - ^ Xt - (vn) 



+ sup 



^r Xt " (vn) - M 



Combining these two inequalities: 



sup \vn (t n ,X l tn ,i) - v n (t n ,X l tn ,i)\ < ^ 
»e4" jeh 



(vn) - ^r Xt " (vn) 



+ sup 

je4 n 



(vn)-^' Xtn (vn) 



Hence, using the triangular inequality, Corollary 3.1, equation (3.30), and the induction hypothesis: 

l + C(T,e) 



sup | v n (t n ,X tn ,i)- v n (t„ ,X tn ,i)\ 



< E n :— C p e 



-pt n 



Mp(T,S,ey ^ 



, r r - P t J + C(T,e) 
+CpC Mp(T,S,s) + En+1 



for some constant C p > which depends only on p. Consequently: 

_ ptn \ + C(T,e) A 1 



where C p > depends only on p. 



1 



hVMp(T,S,ey w 



Mp (T, 5, e) " v2 



□ 



Finally, the combination of Propositions 3.1 3.2, 3.3, 3.5 and 3.6 at time t ~ t proves Theorem 3.1. 
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4 Complexity analysis and memory reduction 



4.1 Complexity 

4.1.1 Computational complexity 

The number of operations required by the algorithm described below is in 0(q 2 .N .M) , where we 
recall that q is the number of possible switches, N is the number of time steps and M is the number 
of Monte Carlo trajectories. 

• The q 2 term stems from the fact that for every i G I g , one has to compute a maximum on j G I q 
(see equation (3.16)). However, this q 2 can be reduced to q as soon as the two following conditions 
are satisfied: 

1. (Irreversibility) The controlled variable can only be increased (or, symmetrically, can only be 
decreased) 

2. (Cost Separability) There exists two functions fci and &2 such that V(t,i,j) € R + x I, x I q , 
k (t,i,j) = k\ (t,i) + hi (t,j). For instance, this is true of afiine costs. 

Indeed, under those two conditions, equation (3.16) becomes: 

vn (t n ,x,i)+ki {t n ,i) = max I hf (t n ,x, j) - fc 2 (£„, j) +E v n (t n+1 ,Xl"f',j ) \ , n = N-l, . . . ,0 

j'Glq , j>i V. L V / J J 



These maxima can be computed in O(q) instead of C(g 2 ) by starting from the biggest element 
i = i q down to the smallest clement i = i\ (in lexicographical order) and keeping track of the partial 
maxima. 

Note that these two conditions hold for the numerical application from Section 5, providing the 
improved complexity O(q.N.M). 



• The N term comes from the backward time induction. 



• The M term corresponds to the cost of a regression, which, in the case of a local basis, can be 
brought down to O (M) (cf. [13]). 



4.1.2 Memory complexity 

The memory size required for solving optimal switching problems (as well as the simpler American 
option problems and the more complex BSDE problems) by Monte Carlo methods is often said to 
be in O(N.M), because, as the Euler scheme is a forward scheme and the dynamic programming 
principle is a backward scheme, the storage of the Monte Carlo trajectories seems inescapable. This 
fact is the major limitation of such methods, as acknowledged in [16] for instance. 

Since such a complexity would be unbearable in high dimension, we describe below a general memory 
reduction method to obtain a much more amenable 0(N + M) complexity (or, more precisely, of 
0(m.N + q.M) with m <§C M). This improvement really opens the door to the use of Monte Carlo 
methods for American options, optimal switching and BSDEs on high-dimensional practical appli- 
cations. Note that this tool can be combined with all the existing Monte Carlo backward methods 
which (seem to) require the storage of all the trajectories. 

A drawback of this tool is that it is limited to Markovian processes. However, one can usually 
circumvent this restriction by increasing the dimension of the state variable. 
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4.2 General memory reduction method 
4.2.1 Description 

The memory reduction method for Monte Carlo pricing of American options was pioneered by [18] 
for the geometric Brownian motion, and was subsequently extended to multi-dimensional geometric 
Brownian motions ([19]) as well as exponential Levy processes ([20]). These papers take advantage 
of the additivity property of the processes considered. However, as briefly hinted in [59], the memory 
reduction trick can be extended to more general processes. In particular, it can be combined with 
any discretization scheme, for instance the Euler scheme or Milstein scheme, as long as the value of 
the stochastic process at one time step can be expressed via its value at the subsequent time step. 

From a practical point of view, the production of "random" sequences usually involves wisely chosen 
deterministic sequences, with statistical properties as close as possible to true randomness (cf. [39] 
for instance for an overview). These sequences can usually be set using a seed, i.e. a (possibly 
multidimensional) fixed value aimed at initializing the algorithm which produces the sequence: 

rand() rand() rand() rand() 
{setseeds} — > e\ — > e-i — > ■■■ — > e n (4.1) 

where the randQ operation consists in going to the next element of the sequence. Now two useful 
aspects can be stressed. The first is that one can usually recover the current seed at any stage of 
the sequence. The second is that, if the seed is set later to, say, once again the seed s from equation 
(4.1), then the following elements of the sequence will be once again E\, £2, ••• In other words, 
one can recover any previously produced subsequence of the sequence (£n) n>1 , provided one stored 
beforehand the seed at the beginning of the subsequence. This feature is at the core of the memory 
reduction method, which we are going to discuss below in a general setting. 

Consider a Markovian stochastic process (X t ) t>0 , for instance the solution of the stochastic differen- 
tial equation (2.2), recalled below: 

X = x e R d 
dX s = b(s,X s )ds + a(s,X s )dW s 

The application of the Euler scheme to this equation can be denoted as follows: 

4 +1 = / (•<•;• ;) (4-2) 

f(x,e) := x + b (ti, x) h + a (tii x) evh (4-3) 

where Vi S [0,iV — 1] and Vj <E [1,-M], e\ € M. d is drawn from a d-dimensional Gaussian random 
variable. Suppose that for any e <E R d , the function x h-> f(x,e) is invertible (call f luv its inverse). 
Then, starting from the final value x\ of the sequence (4.2), one can recover the whole trajectory 
of X : 

4 A, (•<•;....-;) (4.4) 

as long as one can recover the previous draws • ■ ., £q- The following pseudo-code describes an 

easy way to do it. 
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Algorithm 1 Euler Scheme 



Inverse Euler Scheme 



1 % Initialization 

2 for j from 1 to M 

3 X[j] <- xj 

4 end for 

5 

6 % LOOP 1: Euler scheme 

7 for i from to N-l 

8 S [ i ] <- getseed () 

9 for j from 1 to M 

10 E <- rand(d) 

11 X[j] <- f(X[j],E) 

12 end for 

13 end for 

14 S[N] <- getseed () 



1 % LOOP 2: Inverse Euler scheme 

2 for i from 1 down to 

3 setseed (S [ i ] ) 

4 for j from 1 to M 

5 E <- rand (d ) 

6 X[j] <- finv(X[j] ,E) 

7 end for 

8 end for 

9 setseed (S [N] ) 



The first column of Algorithm 1 corresponds to the Euler scheme, with the addition of the storage 
of the seeds. At the end of the first colum, the vector X contains the last values X J T , j — 1, . . . , M. 
From this point, one can recover the previous values Xf, , i = N — 1, . . . , 0, j = 1, . . . , M as done in 
the second column. 

Inside this last loop, one can perform the estimation of the conditional expectations required by the 
resolution algorithm of our stochastic control problem (equation (2.10)). Compared to the standard 
storage of the full trajectories Xf , i = 0, . . . , N, j = 1, . . . , M, the pros and cons are the following: 

• The number of calls to the rand () function is doubled. 

• The memory needed is brought down from O (M x N) to O (M + N) (storage of the vector space 
and the seeds). 

In other words, at the price of doubling the computation time, one can bring down the required 
memory storage by the factor min [M, N) , which is a very significant saving. Moreover, the theoretical 
additional computation time can be insignificant in practice, as the availability of much more physical 
memory makes the resort to the slower virtual memory much less likely. 

Remark 4.1. Even though the storage of the seeds does take O (N) in memory size, the constant may 
be much greater than 1. For instance, on Matlab®, a seed from the Combined Multiple Recursive 
algorithm (refer for instance to [39] for a description of several random generators) is made of 12 
uint32 (32-bit unsigned integer), a seed from the Multiplicative Lagged Fibonacci algorithm is made 
of 130 uint64, and a seed from the popular Mersenne Twister algorithm is made of 625 uint32. 

In order to relieve the storage of the seeds, we now provide a finer memory reduction algorithm 
(Algorithm 2). Although Algorithm 2 requires three main loops, it enables to perform the last loop 
without fiddling the seed of the random generator, and without any vector of seeds locked in memory, 
which will thus be fully dedicated to the regressions and other resolution operations. Moreover, the 
first two main loops can be performed beforehand once and for all, storing only the last values of the 
vector X as well as the first seed S [0] . Finally, if the random generator is able to leapfrop a given 
number of steps, the first loop can be drastically reduced. 
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Algorithm 2 General Memory Reduction Method 



1 


% LOOP 1: Seeds storage 


1 


% LOOP 2: Euler scheme 




2 


for i from to N-1 


2 


for i from to N-1 




3 


S [ i ] <— g e t s e e d ( ) 


3 


setseed (S [N— i -1]) 




4 


for j from 1 to M 


4 


for j from 1 to M 




5 


E <- rand (d ) 


5 


E <- rand (d ) 




6 


end for 


G 


X[j ] <- f (X[ j ] ,E) 




7 

8 


end for 


7 

8 


end for 
end for 




9 


% Initialization 


9 


setseed (S [0] ) ; free (S) 




10 


for j from 1 to M 


10 






11 


X[j] <- xj 


11 


% LOOP 3: Inverse Euler 


scheme 


12 


end for 


12 


for i from N— 1 down to 





13 


% 


13 


for j from 1 to M 




14 


% 


14 


E <- rand (d ) 




15 


% 


15 


X[j] <- finv(X[j] 




16 


0/ 
70 


1G 


end for 




17 


"/ 
70 


17 


end for 





4.2.2 Numerical stability 

Theoretically, the trajectories produced by the Euler scheme (4.2) and the inverse Euler scheme (4.4) 
are exactly the same. In practice however, a discrepancy may appear, the cause of which is discussed 
below. 

On a computer, not all real numbers can be reproduced. Indeed, they must be stored on a finite 
number of bits, using a predefined format (usually the IEEE Standard for Floating-Point Arithmetic 
(IEEE 754)). In particular, there exists an incompressible distance e > between two different 
numbers stored. This causes rounding errors when performing operations on real numbers. 

For instance, consider x € K and an invertible function / : M H> R. Compute y = f (x) and then 
compute x = /; nv (y). One would expect that x = x, but in practice, because of rounding effects, one 
may get x = x + ez for a small e > 0, where z is a discrete variable, which can be deemed random, 
taking values around zero. This phenomenon is illustrated on Figure 4.1, which displays a histogram 
of x — x for n = 10 7 different values of x £ [0, 1] and for the simple linear function / (x) = 2x + 3. 

1 00% r 

80% - 

60%- 

40% - 

20% - 

0%- 
-6 



We now describe how this affects our memory reduction method. Recall equation 4.2: 




-4-2 2 4 6 

-..-.-16 

x 10 

Figure 4.1: Histogram of rounding errors 
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Now, instead of equation (4.4), the inverse Euler scheme will provide something like: 

j j 

Ut N — x t N 

Vi = f™(?i i+1 A)+e4 (4-5) 

for a small e > 0, where zf, i — 0, . . . , N, j — 1, . . . , M, can be deemed realizations of a discrete 
random variable Z, independent of W. The distribution of Z is unknown, but data suggests it may 
be innocuously assumed centered, symmetric, and with finite moments. 

We are now interested in studying the compound rounding error y ti — Xt t as a function of e. Of 
course, its behaviour depends on the choice of / (equation (4.3)). Below, we explicit this error on 
two simple examples: an arithmetic Brownian motion and an Ornstein-Uhlenbeck process. These 
two examples illustrate how the compound rounding error can vary dramatically w.r.t. /. 

First example: arithmetic Brownian motion Consider first the case of an arithmetic Brownian 
motion with drift parameter /j and volatility parameter a . Here / and its inverse are given by: 

f(x,e) = x + nh + aVhs 
/i„v (x, e) = x — fih — aVhe 

Hence, using equation (4.5), for every j — 1, . . . , M: 

N-l 
k—i 

In other words, the compound rounding error behaves as a random walk, multiplied by the small 
parameter e. Hence, as long ast < ft (which is always the case as real numbers smaller than e cannot 
be handled properly on a computer), this numerical error is harmless. 

Remark that a similar numerical error arises from the algorithms proposed in [18] , [19] and [20], 
but, fortunately, as discussed above, this error is utterly negligible. 

Second example: Ornstein-Uhlenbeck process Now, consider the case of an Ornstein-Uhlenbeck 
process with mean reversion a > 0, long-term mean /i and volatility a . Here: 

f(x,e) = x + a (/i — x) h + avhe 
/ itlv (x, e) = - — — - (x - afih - aVhe^j 

Using equation (4.5), for every j = 1, . . . , M the compound error is given by: 

vi - x t 

As (1 — ah) N ~ exp (qT) when h — > 0, one can see that, as soon as T > — , this error may 
become overwhelming. This phenomenon is illustrated on Figure 4.2a on a sample of 100 trajectories. 

In order to mitigate this effect, we propose to modify the Algorithm 2 as follows: in its second loop 
(usual Euler scheme), instead of saving only the last values x J Tl one may define a small subset IIcII 
and save the intermediate values x J t ., ti £ 5. Then, in the last loop (inverse Euler scheme), every 
time that U £ H, the current value of the set x\. may be recovered from this previous storage. 

Figure 4.2b illustrates the new behaviour of the compound rounding error with this mended algorithm, 
on an example with T = 10 years and 4 intermediate saves (in addition to the final values). 

The drawback of this modification, of course, is that it multiplies the required storage space by the 
factor #11. However, this remains much smaller than the O (M x N) required by the naive full 
storage algorithm. 



JV-l 



ahf 
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5 Application to investment in electricity generation 

This section is devoted to an application of the resolution method studied in Section 2. We choose to 
apply it to an investment problem in electricity generation on a single geographical zone. We intend to 
show that it is possible to provide a probabilistic outlook of future electricity generation mixes instead 
of a deterministic outlook provided by planification methods. Nevertheless, the problem presents so 
many difficulties that addressing all of them in the same model is unresonable. Some aspects have 
thus to be left aside. Our goal here is to show that the algorithm described in Section 3 can handle 
high-dimensional investment problems. We focus on the influence of investment decisions on the spot 
price, consistently with the fundamentals of the electricity spot price formation mechanism. 

Although the strategic aspect of investment is an important driver of utilities' decisions, this aspect 
is beyond the scope of our modeling approach. There exist models limited to a two-stage decision 
making (see for instance [48]), but in the case of continuously repeated multiplayer game models, 
defining what is a closed-loop strategy is already a difficulty (see Sec. 2 in [7]). 

We did not consider time-to-build in this implementation either. Relying on the fact that it is possible 
to transform an investment model with time-to-build into a model without time-to-build by replacing 
capacities with committed capacities (see [9, 1] for implementations in dimension one, and [26] in 
dimension two), we left this aspect for future work. 

Finally, we did not consider the dynamic constraints of power generation. Their effect on spot prices 
is well-known: they tend to increase spot prices during peak hours and to decrease them during 
off-peak hours (see [41]). However, we assume here that this effect is negligible compared to the 
effect induced by a lack or an excess of capacity. 

Thus, we focused on the following key factors of electricity spot prices: demand, capacities (including 
random outages) and fuel prices. Our model is based on [5, 4], where the electricity spot price is 
defined as a linear combination of fuel prices multiplied by a scarcity factor. This model exhibits the 
main feature wanted here, which is that the spot price, being determined both by the fuel prices and 
the residual capacity, is directly affected by the evolution of the installed capacity. When the residual 
capacity tends to decrease, spot prices will tend to increase, making investment valuable. Thus, in 
this model, investments are undertaken not on the specific purpose of satisfying the demand but as 
soon as they are profitable. In our example, new capacities are invested according to the criterion 
of value maximization. Energy non-served and loss of load probability may still be adjusted through 
the price cap on the spot market. 

In this section, we first detail the chosen modelling and objective function (which will be shown to be 
encompassed in the general optimal multiple switching problem (2.1)), and then solve it numerically 
using the general algorithm developed in the previous sections. 
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5.1 Modelling 



The key variable in order to describe our electricity generation investment problem is the price of 
electricity. More precisely, the key quantities are the spreads between the prices of electricity and 
other energies. To model these spreads accurately, it may be worth considering a structural model for 
electricity (cf. the survey [14]). Here we choose such a model, mainly inspired by those introduced in 
[5] and [4], albeit amended and customized for a long-term time horizon. All the variables involved 
are detailed below. 



5.1.1 Electricity demand 

The electricity demand, or electricity load, at time t on the given geographical zone considered is 
modelled by an exogenous stochastic process (D t ) t>0 : 

D t = f {t)+Z° t (5.1) 

where Z° is an Ornstein-Uhlcnbeck (henceforth O.U.) process: 

dZ° = -a Z°dt + P dW t D 

where ao and /3o are constants, and fo is a deterministic function that takes into account demand 
seasonalities: 

fo (t) = d l + d 2 cos ^r^j-^ + f week (t) (5.2) 

where dj, 1 < j < 3 are constants, and, assuming that t is expressed in years, l\ = 1 (yearly 
seasonality) , and f wee k is a periodic non-parametric deterministic function describing the intra- week 
load pattern. 



5.1.2 Production capacities 

Let d! be the number of different production technologies. Denote as It = (if, ...,1^ the installed 
production capacities at time t. They represent the maximum amount of electricity that is physically 
possible to produce. These fleets can be modified: at a given time t„, one can decide to build (or 
dismantle) an amount £ n of capacities: 

Ir n = I T - + Cn , n > (5.3) 

Denote as a = (r„, C«)„>i the corresponding impulse control strategy, where (T n ) n>0 is an increasing 
sequence of stopping times with r„ /• oo when n — > oo, and (Cn)„>o i s a sequence of vectors 
corresponding to the increases (or decreases) in capacities. Apart from these variations, /( will be 
deemed constant, i.e.: 

I t =I -+ ]T C». (5.4) 

n , r n <t 

Now, denote as C t = (cj, . . . ,Cf^j the available production capacities. Because of spinning reserves, 
maintenance and random outages, these quantities are lower than the installed capacities It, which 
represent their physical maximum. In other terms, C't is a fraction of It'- 

Ci = P t x4 (5.5) 

for every 1 < i < d 1 , where A\ corresponds to the rate of availability of the i th production technology. 
Therefore one must choose a model for the process A t that ensures that it stays within the interval 
[0,1]. 

One possibility would be to model it as a Jacobi process (see for instance [58] ,where it is used to model 
stochastic correlations, and the references therein for more information on this process). This process 
is however tricky to estimate and simulate (see [31] for the description of some possible methods), 
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and its main simulation method (the truncated Euler scheme) disables our memory reduction method 
described in Subsection 4.2. Hence we look for a simpler model. 

In [60], a detailed structural model for electricity is developed, which includes renewable energies like 
wind and solar. In particular, wind power infeed efficiency (which belongs to [0, 1]) is modelled as 
a logit transform of an Ornstein-Uhlenbeck process with seasonality. Adapting this idea, we model 

( A l)l>l~ d as follows: 

Al:=T(f t (t) + Zl) (5.6) 
where Z, f and T are chosen as follows: 

• Z % is an O.U. process : 

dZ\ = -aiZidt + PidWf 
where 014 > 0, fii > and ( W? ) is a Brownian motion. 

• The deterministic function fi accounts for the seasonality in the availability of production capac- 
ities: 

/ i (t)=4+4cos^27r^-^^ (5.7) 

where c\, 1 < k < 3, 1 < i < d' are constants. This seasonality stems from the maintenance 
plannings, which usually mimic the long term seasonality of demand, which in turn originates in 
the seasonality of temperature. 

• The function T ■ K — > [0,1] is here to ensure that Vf > 0, A t € [0, 1] . One can choose the 
versatile logit function as in [60], or any other mapping of K into [0, 1]. For instance, any cumu- 
lative distribution function would be suitable. As the process Z is Gaussian and asymptotically 
stationary, we choose for T the (standard) normal cumulative distribution function, as it makes, 
in particular, the calibration process trivial. 



5.1.3 Fuels and CO2 prices 

For each technology i, denote as S} the price of the fuel i to produce electricity at time t. In the 
particular case of renewable energies, which, per se, do not involve traded fuels, the corresponding 
S} can be chosen to be zero. Moreover, define S® as the price of C0 2 - Denote as S t the full vector 

[St > &t i • ■ • i &t ) ■ 

Now, we introduce the multiplicative constants needed to convert theses quantities into €/MWh. 
For each technology i — 1, . . . , d 1 , let hi denote its heat rate, and h® denote its CO2 emission rate. 
Hence, the quantity 

§1 := h?X + hiSt (5.8) 

expressed in €/MWh, corresponds to the price in € to pay in order to produce lMWh of electricity 
using the ith technology. We note h° = (h%,...,h%) G R d ' and h = (hi, . . . , h d >) £ R d ' . 

Remark 5.1. One can choose to add a fixed cost into the definition of S\. This is all the more so 
relevant for technologies whose fixed costs outweigh the cost of fuel (e.g. nuclear). 

Adapting the work of [10], we model St as a multidimensional, cointegrated geometric Brownian 
motion: 

dS t = ES t dt + EStdWf 

where S and X are {d 1 + 1) x {d 1 + 1) matrices with 1 < rank (S) < d', and (W / t S ) t>0 is a (d' + 1)- 
dimensional Brownian motion. This model ensures the positivity of prices, as well as the existence 
of long-term relationships between energy prices (the relevance of which is illustrated, for instance, 
in [49]). 
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5.1.4 Electricity price 



We model the price of electricity using a long-term structural model. First, we define the marginal cost 
of electricity using the previously introduced variables. For any time t > 0, define the permutation 
(1) , . . . (M) of the numbers l,...,M, such that sf 3 < . . . < S { t M) . Then, define cf } as the total 

capacity available at time t from the i first technologies, i.e. cf := J2j<i • Using these notations 
and equation (5.8), the marginal cost of electricity at time t is given by: 

M-i 

MC t : = S«l {A < cf '} + £ Sf 1 {C [ r 1] < A < cf) + sl M h < D t ) 

i=2 

M-l 

= + £ - § t (i) ) i {a - > o} 

i=l 

Refer to [5] for more details on marginal costs. Remark that the price of CO2 emissions is explicitly 
included in the marginal cost (through equation (5.8)). 

Now, we are going to use this marginal cost as a building block of our price model, along with some 
power law scarcity premiums (along the lines of [4]) as well as a fixed upper bound . 

First, consider two points (£1,1/1) and (x2,y 2 ) in R 2 - One can always find three positive constants 
a := a (xi, X2, Hi, 2/2)1 b := b (x\, x<i , 2/1, 1/2) and c :— c (xi, X2, 2/1, 2/2) such that the function: 

p{x) :=p(x;xi, x 2 , 2/1,2/2) = 7— he (5.9) 

— x 

satisfies p(x\) = 2/1 and p{x2) — 2/2 2 - 

Using this notation, introduce the price Pt of electricity, defined as follows: 

p t := 5«i{A<0} + {^+p(A;0,cf\sW 5f))}i{o<A<cf ) } 



e +p (A^r^u^p 15 ) } 1 {^r 1 ' < a < c«} 

2 = 2 




(5.10) 



where M max > is a fixed upper bound on the price of electricity. In particular, the last term, the 
one involving M max , enables price spikes to occur (when the residual capacity is small). 

Moreover, thanks to the knitting function (5.9), the electricity price P is a Lipschitz continuous 
function of the structural variables D, C and S 3 , which is what motivated this specific choice of 
model. 



5.1.5 Objective function 

We now explicit the objective function of the investor in electricity generation. Suppose that, at time 
t, an agent (a producer, or an investor) modifies the level of installed capacity of type j £ 
from If_ to I 3 S = P{_ + £ 3 , s > t . It generates the cost: 

r„/+ + c^f .c -o 

k (C J ) := I , C J = 

1 Indeed, in the French, German and Austrian markets for instance, power prices cannot be set outside the 
[—3000, 3000]€/MWh range, see http://www.epexspot.com/en/product-info/auction.. 

2 For instance, fix a > 0, then define b = ~ Aei + £2 + »J (x% — £'l) 2 + 4a and finally c = yi — gd|rp 

3 Rigorously, this property requires that C does not reach zero. One can, for instance, add a fixed minimum 
availability rate 1 2> a m in > to the definition (5.6), replacing T by a m i n + (1 — a m i n )T 
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where and are the fixed and proportional costs of building new plants of type j, and 
and kP~ are the fixed and proportional costs of dismantling old plants of type j. 

Consider the case of new plants (( 3 > 0). Assuming that the global availability rate(5.6) of technology 
j applies to the new plants, they can then produce up to ( 3 A 3 S , s > t, or, more precisely, according 
to the stack order principle: 



assuming that, in the stack order, the new plants are called before the older plants /(_ of the same 
technology (as they can be expected to have an at least slightly better efficiency rate compared to 
the older plants of the same technology, a phenomenon that can be seen as partly captured by the 
function (5.9)). 

At time s > t, this production is sold at price P s , but costs S s to produce (if P s < S s , then of course 
the producer chooses not to produce). In addition, regardless of the output level, there may exist a 
fixed maintenance cost Kj . Summing up all these gains, discounted to time t using a constant interest 
rate p > 0, the new plant yield a revenue of: 



e'" s (min {c'Aj, D s - C ( s j 1] } x (p s - - n t 



ds 



(noticing that with our power price model, j.D s — C^J ^ < o| |p s — S 3 < oj). This was the 

cost-benefit analysis for one quantity Q 3 of new plants. Now, consider as a whole the full fleet of the 
geographical zone considered. Maximizing the expected gains along the potential new plants yields 
the following value function: 



v (t, x, i) = sup E 

aeAt.i 



d pOG / _|_ \ 

J2 j t e ~ ps ( min { c i> D > - } x { p s - si) - «ij ds - ^ pTnk 



(c 



(5-11) 

where the strategics a affect the installed capacities (equations (5.4)), hence also the available ca- 
pacities (equation (5.5)) as well as the power price (equation (5.10)), and where the cash flows are 
purposely discounted up to time 0, the time of interest. 

Remark 5.2. Replacing P in (5.11) by its definition (5.10), it is patent that this objective function 
fits into the mould studied thoroughly in Section 3. In Subsection 5.2 below, our algorithm will be 
applied to this specific objective function. 

Remark 5.3. Remark that under this modelling, the demand is satisfied as long as it does not exceed 
the total available capacity. Indeed, the effective output of the plant Q 3 is equal to 1 |p s — S 3 S > j x 

min |<^M^, (^D s — ^ |. It is indeed governed by the electricity spot price level, but, as under 

our modelling 1 |p s — S 3 > oj = 1 |l) s — C^J ^ > oj, summing up the effective outputs of all the 

power plants yield £)f = i min {cj, D s - C^~ 1] } x 1 {d s -C { ^ 1] > o} = min |d s ,C^'^|. 



5.2 Numerical results 



Finally, we solve the control problem described in Subsection 5.1 on a numerical example, using the 
algorithm detailed in Subsection 3 combined with the general memory reduction method described 
in Subsection 4.2. 

Our purpose here is not to perform a full study of investments in electricity markets, but a more 
modest attempt at illustrating the practical feasibility of our approach, with some possible outputs 
that the algorithm can provide. 
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We consider a numerical example including two cointegrated fuels (in addition to the price of CO2): 
one "base fuel" and one "peak fuel", starting respectively from 40€/MWh and 80€/MWh. Hence, 
using the notations from Subsection 5.1, d! = 2 (two technologies) and d = 6 ( electricity demand, 
CO2 price, two fuel prices and two availability rates). The main choices of parameters for this 
application (initial fuel prices and volatilities, initial fleet and proportional costs of new power plants) 
are summed up in Table 5.1. Moreover, the demand process starts from Dq — 70GW and does not 
integrate any linear trend. 



i 


oi 


a 1 




< + 


1 


40€/MWh 


5% 


67GW 


0.24 10 y €/GW 


2 


80€/MWh 


15% 


33GW 


2.00 10 9 €/GW 



Table 5.1: Model parameters 



In order to take into account the minimum size of one power plant we restrict the values of the 
installed capacity process(5.4) to a (bi-dimensional) fixed grid A d , with a mesh of 1GW. We make 
the simplifying assumptions that investments are irreversible, and that no dismantling can occur 
(recall from Subsection 4. 1 the computational gain provided by this assumption) . 

Remark 5.4. If such a grid is indeed manageable in dimension d' = 2, it may less be the case if 
additional technologies were considered. However, as discussed in [55] equation (3.2), instead of 
performing one regression for each i S A d , one can solve equation (3.16) at time i, by only one 
(d + d' )-dimensional regression, by choosing an a priori law for the randomized control £ t .. The error 
analysis from Section 2 can be easily generalized to such regressions in higher dimension. 

Finally, we consider the following numerical parameters. We choose a time horizon T = 40 years and 
a time step h = (i.e. two time steps per day, allowing for some intraday pattern in the demand 
process) but allow for only one investment decision per year. For the regression, we consider a basis 
of b = 2 d — 64 adaptative local functions, chosen piecewise linear on each hypercube (which is a 
bit more refined than the piecewise constant basis studied in Section 3) on a sample of M = 5000 
trajectories. 

The numerical results obtained under this set of parameters are displayed on Figures 5.1 and 5.2. 

First, Figure 5.1 deals with the optimal strategies. Figure 5.1a displays the time evolution of the 
average as well as the variability of the optimal fleet (only the new plants are shown). One can 
distinguish a first short phase characterised by the construction of several GW of peak load assets, 
followed by a much slower second phase involving the construction of both base load and peak load 
assets. Moreover, the variability of the optimal fleet increases over time. The detailed histogram of 
the optimal strategy at time T = 40 years is displayed on Figure 5.1b, where it is combined with 
the price of fuel. One can see that the more the peak fuel is expensive (and hence both fuels are 
expensive on average, as they are cointegrated), the more constructions of base load plants occur. 
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The fact that the average fleet seem to converge is related to the fact that this numerical example 
does not consider any growth trend in the electricity demand (see equation (5.2)). Otherwise, more 
investments would occur, indeed, to keep the pace with consumption. 

Then Figure 5.2 provides information on the price of electricity. Figure 5.2a displays the time 
evolution of the electricity spot price density. For better readability, each density covers one whole 
year. One can see how the density moves away from the initial bimodal density (with prices clustering 
around the initial prices of the two fuels) towards a more diffuse density. Moreover, the downward 
effect of investments on prices can be noticed. This downward effect is even more visible on Figure 
5.2b. It compares the effect on electricity prices of three different strategies: the optimal strategy, 
the optimal deterministic strategy (computed as the average of the optimal strategy) , and the do- 
nothing strategy. For each strategy, the joint time-evolution of the yearly median price and the yearly 
interquartile range are drawn. As expected, prices tend to be higher and more scattered without 
any new plant. Nevertheless, on this specific example, the price distribution under the optimal 
deterministic strategy is close to that under the optimal strategy (only slightly more scattered) . 




Median Price (Euro/MWh) 

(a) Time evolution of electricity spot price density (b) Comparison between investment strategics 

Figure 5.2: Electricity spot price 



These few pictures illustrate the kind on information that can be be extracted from the resolution of 
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this control problem. Of course, as a by-product of the resolution, much more can be extracted and 
analyzed (distribution of income, CO2 emissions, optimal exercise frontiers, etc) if needed. 



6 Conclusion 

In this paper, we presented a probabilistic method to solve optimal multiple switching problems. We 
showed on a realistic investment model for electricity generation that it can efficiently provide insight 
into the distribution of future generation mixes and electricity spot prices. We intend to develop 
this work in several directions in the future. First, we wish to take into account more generation 
technologies, most notably wind farms, nuclear production, as well as solar distributed production. 
These additions would raise the dimension of the problem from eight to fifteen. Yet another range of 
innovations in numerical methods will be necessary to overcome this increase in dimension. Second, 
we wish to take time-to-build into account. And last but not least, we wish to adapt the problem to 
a continuous-time multiplayer game and contribute to the quest for an efficient algorithm to solve it. 
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